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We have performed molecular dynamics simulations of solid-solid contacts lubricated by a model 
fluid displaying many of the properties of water, particularly its expansive freezing. Near the region 
where expansive freezing occurs, the lubricating film remains fluid, and the friction force decreases 
linearly as the shear velocity is reduced. No sign of stick-slip motion is observed even at the lowest 
velocities. We give a simple interpretation of these results, and suggest that in general good boundary 
lubrication properties will be found in the family of materials with expansive freezing. 



The understanding from a microscopic point of view of 
the tribological properties of materials in the boundary 
lubrication regime is of great importance for the tech- 
nological development of new lubricants for new appli- 
cations Micro-mechanical devices and high density 
magnetic storage media are just two examples where the 
detailed properties of the lubricants may affect in an 
important way the performance of the system. To at- 
tain low friction coefficients with lubricating films a few 
molecular layers thick is just the most obvious goal of 
lubricant design. 

Fluid lubrication is a widespread strategy to reduce 
friction. However it suffers of some limitations in the 
boundary regime |^ , where the fluid is only a few molecu- 
lar layers thick. Fluids have a general tendency to attain 
a solid-like structure when squeezed between two solid 
surfaces [^||^. This confinement-induced solidification is 
responsible for microscopic stick-slip motion (SSM) |^] 
when the confining surfaces are sheared past each other. 
SSM in turn generates a finite friction force in the limit 
of low velocities, and then a non-zero friction coefficient. 
In fact, mechanisms that prevent SSM have been seen to 
reduce dramatically the friction force between the lubri- 
cated surfaces . 

In the boundary regime a lubricant composed of com- 
plex, branch chain molecules usually performs better 
than one formed by spherical, or single chain molecules, 
as the former does not readily arrange in a solid-like 
structure Q, and allows for a smooth shearing of the 
surfaces, in a state of very low friction. In addition, 
fiuids with branch chain molecules have typically higher 
bulk viscosity than those with spherical, or single chain 
molecules. This justifies the experimentally observed (see 
Table 1 in Ref. |^ ) negative correlation existing between 
bulk viscosity and friction coefficient under boundary lu- 
brication conditions. 

Water is a noteworthy exception to this rule The 
water molecule is a small, almost spherical particle, and 
the viscosity of water is very low. However, carefully 
controlled experiments using the surface force apparatus 
show that water is an excellent boundary lubricant at a 
microscopic scale (7[j^. In fact, the friction coefficient be- 



tween two water lubricated, atomically flat mica surfaces 
was found to be around 0.01, in the range of the 
lowest values obtained with other fluids in the boundary 
regime. The reason for the anomalously good lubrication 
properties of water has not been addressed in detail. In 
this paper we show by studying a model system, that 
the good lubrication properties of water can be associ- 
ated to its expansive freezing, and that good boundary 
lubricants can be fetched in general within the class of 
materials having this property. 

Many of the anomalous properties of water (and tetra- 
hedrally coordinated liquids in general [|o|) have been 
recently systematized by the use of a model system of 
spherical particles interacting through a soft-core poten- 
tial [|l]-|l|. Its main characteristic is the existence of two 
stable relative distances between particles, as a function 
of the applied pressure. It has been suggested that these 
two distances are realized in water owing to two distinct 
arrangements of hydrogen bonds |T^. We use one such 
model jm to simulate a lubricant film. The interaction 
potential as a function of the interparticle distance r has 
a strict hard core at a distance r = cr and a repulsive 
ramp extending from a to aa {a > I, a — 1.5 hereafter), 
in which the potential energy changes linearly from e to 
zero. Natural units for this problem are the length a, the 
energy e, and the mass of the particles m, from which a 
time scale r = ■^/mo^/e is defined. 

The bulk system freezes into a closed packed lattice 
with lattice parameter ~ aa at low temperatures and low 
pressures. The melting curve of this crystalline structure 
was determined previously by Monte Carlo simulations 
|lT| , and is reproduced in Fig. 1(a) for completeness. 
The melting temperature increases as a function of pres- 
sure at the lowest values of P and T, consistent with the 
fact that in this region particles behave as hard spheres 
of radius aa/2. At higher pressures, the melting line 
reenters toward lower temperatures as neighbor particles 
can get closer than aa Jisf. Other crystalline structures 
exist at higher pressures | ]l4P^ - 

The geometrical arrangement we use for the simula- 
tions is similar to that used in Refs. [17] and [18], and is 
sketched in Fig. |l|(b) . It qualitatively models the config- 
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uration used in the surface force apparatus Two op- 
posing pieces of crystalline fee solids exposing (111) faces 
are lubricated by our model system. Particles within the 
solids are assumed to be fixed to their lattice positions, 
and the interactions between particles of the fluid and 
those of the solids are taken equal to the intra-fluid in- 
teractions. The lattice parameter of the solid blocks uq is 
chosen for each pressure as that of the bulk solid at the 
melting temperature. This is not a crucial choice, and is 
made only to avoid introducing a bias between the dif- 
ferent simulated points. The lubricant particles within 
the gap between the solids are at dynamic equilibrium 
with lubricant particles outside the gap. For the static 
(no shear) simulations, periodic boundary conditions are 
applied in the three directions. The solid blocks are of 
finite size Ix in the x-direction, and extend over the whole 
length of the simulation cell ly in the y-direction. The size 
of the simulation box in the a;-direction is substantially 
larger than l^, allowing to reach a dynamical equilibrium 
between the lubricant inside the gap, and the bulk fluid 
outside it. Ix is allowed to change according to the value 
of a hydrostatic external pressure P. Real fluids always 
possess some attraction between particles, that we have 
not included. It is known however that the attraction 
between particles (both fluid-fluid and fluid-solid) in the 
van der Waals limit can be readily included in the results 
obtained without attraction by the renormalization of the 
external pressure P |]ll| . Our hydrostatic pressure P ac- 
counts for the attraction between particles in the van der 
Waals limit, this attraction being of the same strength 
between fluid-fluid particles, and fluid-solid ones. Other 
cases remain to be studied. 

We have chosen x ly equal to 15 x 15cr^. Each lubri- 
cant monolayer within the gap has about 110 particles, 
and the total number of fluid particles simulated is 800. 
We have performed checks with systems of up to 3500 
fluid particles and x ly equal to 30 x SOa^, to assure that 
size effects do not spoil the results. We simulated systems 
at points A and B, indicated in Fig. 1 {Pa — 0.1 e/a^, 
Ta = 0.035 e/fcs, Pb = 0.9 e/cr^, Tb = 0.13 e/fcs). We 
calculated the normal force exerted by the lubricant -the 
solvation force- and studied the structure of the film as a 
function of the distance d between the two solid surfaces. 
In Fig. 2 we show results at point A, where the system 
behaves basically as a hard sphere system. When d is 
only a few molecular radii, the solvation force per unit 
area Fs shows a well known |0,^-22| oscillatory behav- 
ior (centered at the value of the hydrostatic pressure P) 
associated to the layering of the fluid. The number of 
particles inside the gap Up (which is normalized to one 
for a full layer of particles) has a staircase-like structure, 
remaining approximately constant as increases (we 
describe the process as d is reduced), and jumping down 
coincidentally with the reduction of Fs, when a layer of 
lubricant is squeezed out of the gap. Gray regions in Fig. 



2 indicate precisely these squeezing-out regions, where 
the number of layers is changing. We measure the hor- 
izontal order within the lubricant through the structure 
factor / of the layers adjacent to the solids (calculated at 
\k\ — 47r/ao), and the two dimensional (horizontal) diffu- 
sivity of particles in the gap D. Within a plateau of np, 
/ increases and D decreases as the film is compressed, 
indicating that the film becomes more structured as F^ 
increases. In the squeezing-out regions (the gray regions) 
/ drops and D maximizes, indicating that the film be- 
comes less structured and more fluid. This whole picture 
is consistent with previous studies 18|, and it can also 



be given a straightforward interpretation: if we assimi- 
late Fs to a sort of local pressure in the film, the increas- 
ing of order as increases is consistent with the fact 
that for the bulk material positional order increases and 
diffusivity decreases with the applied pressure. When a 
layer is being squeezed out, more space is available to the 
remaining particles that then attain a more fluid conflg- 
uration. 

Results are shown in Fig. 3 near the region where ex- 
pansive freezing occurs (point B in Fig. 1). Fg shows the 
same oscillating behavior as before. However, the struc- 
ture of the fllm is now qualitatively different. The film 
remains in a liquid-like state as F^ increases (as seen from 
the values of / and D). This occurs because particles, 
being at distance 1.5cr from their neighbors, can go 
closer (up to distance a) upon a reduction of d, without 
generating a higher degree of order in the system. Only 
after Fs decreases the film becomes more solid-like (with 
higher / and lower D). This is also consistent with the 
fact that this kind of model (as well as real water [p3[) 
is known to have a range of pressure in which diffusivity 
increases as pressure increases Close to this region 
(which occurs near the maximum of the P{T) melting 
line), a positive correlation between the values of Fs and 
D can be expected. 

The inverted alternate at point B of liquid-like and 
solid-like structure of the fllm upon compression repre- 
sents a crucial difference from an experimental point of 
view, since typically a condition of constant load is ap- 
plied to the contact. For stability reasons, this external 
load will be sustained by a lubricating fllm at a point 
with negative derivative of Fs with respect to d, namely, 
it will be a solid-like film at point A, and a liquid-like film 
at point B. As we will see, this difference has profound 
consequences in the lubricating properties of the system. 

Shearing simulations were performed in the same ge- 
ometry than that used for the previous calculations, but 
now the upper solid block (which is now allowed to move 
in the y direction) is given a mass M (M = 450 rri), 
and is attached to a horizontal spring (spring constant 
k = 0.45 s/a^) whose free end is moved in the y direc- 
tion at a constant velocity v (see also [4]). An external 
load per unit area L of 0.05 e/a'^ in excess of the hydro- 
static pressure P within the lubricant is applied to the 
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block along the z direction, and the distance d is allowed 
to change dynamically during the simulation. 

Both at points A and B a two-layer thick film was 
first stabilized (corresponding approximately to points 
labeled 1 in Figs. 2 and 3, we will refer to these condi- 
tions as Ai and Bi respectively). The friction force Ff as 
a function of time was calculated for different shear veloc- 
ities, and the results are shown in the main panel of Fig. 
4. We show the time dependence of Fj for the points of 
lowest velocities in the bottom part of Fig. 4. At Ai the 
mean friction force goes to a constant when v is reduced, 
as SSM sets in, around v ^ 0.05 cr/r, which reproduces 
previous results in related models (2^,^. The limiting 
value of the friction force corresponds in our case to a fric- 
tion coefficient /i ~ 0.2. At Bi the friction force decreases 
linearly as v is reduced. No SSM is observed in this case 
|p6| . This is consistent with the fluidity of the lubricant 
even in this highly confined geometry, a fact recently ob- 
served in water The linear dependence the friction 
force Ff{v) on v can be used to estimate an effective vis- 
cosity r]"^ for the film as rj"^ = Ff{v)d/v ~ 0.7 er/a^, 
that is within three times the bulk value obtained by 
an independent measurement of the bulk single particle 
diffusivity (D^"!"^ ~ 0.022 a^/r) and use of the Stokes- 
Einstein relation. This remarkable fact is to be compared 
with the finding that viscosity of water confined to molec- 
ular thicknesses has been shown Q to be no more than a 
few times the bulk value, and it can also be traced back 
to the high value of diffusivity obtained in the static cal- 
culation at point 1 in Fig. 3, which is close to the bulk 
value. 

Since temperature and pressure are different at points 
A and B, an additional check of the proposed scenario 
is desirable. Therefore, we did additional simulations 
in which the behavior is expected to be reversed: we 
performed fixed d simulations at values d = 4.8 a at A, 
and d = 3.7 cr at B (labeled 2 in Figs 2 and 3, we will 
refer to these conditions as A2 and B2 respectively. Note 
that constant load simulations at these values of d would 
be unstable). We have obtained in fact a big reduction 
of the friction force at A2 with respect to Ai, with no 
SSM detectable even at the lowest velocities (see bottom 
of Fig. 4), and an increase of the friction force at B2 
with respect to Bi, where now SSM becomes observable 
at the lowest velocities. 

Summarizing, we have shown that materials with ex- 
pansive freezing are expected to have good boundary 
lubrication capabilities, essentially due to the fact that 
they do not show a strong tendency to solidify as they 
are squeezed. This finding provides a novel strategy in 
the design of materials for boundary lubrication. At the 
same time our results give a qualitative explanation for 
the good lubrication properties and low viscosity of very 
thin films of water. They also predict that experiments 
performed at constant d, at the squeeze-out points (such 
as B2) will produce larger friction forces than constant 



load experiments, contrary to what is expected in nor- 
mal fluids. Similar results are likely obtained in other 
materials with expansive freezing. 
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From the shear simulations wo cannot rule out the exis- 
tence of SSM al lower velocities, but from the results in 
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cient n at Bi is obtained, namely 20 times lower than 
for Al. 
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FIG. 1. (a) The melting line of the low pressure bulk crys- 
talline structure of our model fluid. The points A and B 
are those simulated, with the results shown in Fig. 2 and 
3 respectively, (b) A sketch of the geometrical configuration 
used in the simulations. Fluid particles fill up the simulation 
box. Crystal structure of solid blocks is not shown here. The 
external hydrostatic pressure P is indicated. 




FIG. 2. Structure of the lubricant film at point A of the 
phase diagram as a function of the gap width d. We show 
the equilibrium values of the solvation force Fs, the number 
of particles in the gap Up (normalized to a full layer of par- 
ticles), the structure factor / of the lubricant layers adjacent 
to the solids, and the horizontal diffusivity D of particles in 
the gap {D is measured relative to the bulk three dimensional 
value at point A, D^""^ ~ 0.012 /t). In the gray regions 
there are drops of Fs , associated with the expulsion of a whole 
lubricant layer, as seen by the evolution of Up. This expulsion 
is correlated with a decrease of / and an increase of D, which 
indicates that the film becomes more disordered, or fluid-like, 
in these regions. 
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FIG. 3. Same as Fig. 2 but at point B of the phase diagram. 
The trends for F^ and Up are the same as those seen in Fig 
2, but now the expulsion of a lubricant layer (within the gray 
regions) produces an ordering of the film, as seen from the 
increase of / and the decrease of D (D|""^ ~ 0.022 a^/r). 
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FIG. 4. Top: friction force per unit area Ff as a function 
of shear velocity for points A and B of the phase diagram. 
Continuous lines are results at a constant load L of 0.05 a^/e 
in excess of the hydrostatic pressure (points labeled 1 in Figs. 
2 and 3). Dotted lines are the results of simulations at con- 
stant d, chosen in such a way that now the film is quasi-liquid 
at Aa, and quasi-solid at B2 (point labeled 2 in Figs. 2 and 
3). Bottom: Ff as a, function of time for the points of low- 
est velocities, displaying SSM at Ai (for which Ff goes to a 
finite value as v is reduced) and uniform motion at Bi (for 
which Ff decreases linearly as v is reduced). SSM is strongly 
suppressed at A2 as compared to Ai, and reappears in B2, 
as compared to Bi. 
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